0. STUDY CONTEXT AND OBJECTIVES

Tile: Telepressure & Suicide in Secundary School Spain - Network Analysis: A Network Analysis Approach

Authors:
José-Manuel Guerra1,
Jesús Martínez-Calvo2,
Diego Diaz-Milanes2,3

Affiliation:

1 Department of Social Psychology, University of Seville, Seville, 41018, Spain

2 Department of Quantitative Methods, Universidad Loyola Andalucía, Sevilla, Spain

3 Health Research Institute, University of Canberra, Canberra, ACT, Australia

Reproducibility Statement

This R Markdown file contains the de-identified data and all R code used in the study, ensuring full transparency and reproducibility.

Analysis Workflow

The workflow comprises ten steps:

  1. Data import and preprocessing
    Load the dataset, screen missingness/outliers, and prepare variables.

  2. Descriptive statistics
    Summarize sociodemographics, inspect the distribution of variables, and explore gender differences where relevant.

  3. Undirected network modeling
    Estimate a partial correlation (GGM) network of WTM items using EBICglasso.

  4. Network invariance testing
    Compare network structure, global strength, and specific edges across gender groups via the Network Comparison Test (FDR-adjusted edgewise tests).

  5. Directed network modeling
    Learn a Bayesian Network (DAG) for potential directional dependencies among WTM items with bootstrap aggregation.

  6. Predictive evaluation
    Assess model predictability (e.g., \(R^2\), RMSE) via cross-validation and compare performance across models.

1. DATA LOADING & PROCESSING

In this section, the packages, function design, and data used are listed. Then, data processing—including variable selection and filtering by units of analysis—was performed.

1.1. Import Packages & Functions

# ---- Load libraries (already attached above; keep here if you knit this chunk standalone)
suppressPackageStartupMessages({
  library(haven); library(tidyr); library(dplyr); library(corrplot); library(psych)
  library(parameters); library(ggplot2); library(semTools)
  library(tibble); library(purrr); library(effectsize); library(EGAnet); library(qgraph)
  library(NetworkComparisonTest); library(bootnet); library(mgm); library(bnlearn); library(caret)
  library(forcats); library(knitr); library(kableExtra)
})

# ---- Helper: safe version + loaded flag
safe_version <- function(pkg) {
  v <- tryCatch(as.character(utils::packageVersion(pkg)), error = function(e) NA_character_)
  if (is.na(v)) "" else v
}
is_loaded <- function(pkg) pkg %in% loadedNamespaces()

# ---- Packages + key functions actually used/reported
pkgs <- tibble::tibble(
  Package   = c(
    "readxl","tidyr","dplyr","psych","semTools","ggplot2","effectsize",
    "EGAnet","qgraph","NetworkComparisonTest","bootnet","mgm","bnlearn",
    "caret","forcats","knitr","kableExtra","lavaan","semPlot","corrplot","haven","parameters"
  ),
  Functions = c(
    "`read_excel`",
    "`pivot_longer`, `pivot_wider`, reshaping support",
    "`select`, `mutate`, `filter`, `group_by`, `summarise`, `bind_rows`",
    "`describe`, `polychoric`, `alpha`, `omega`, `splitHalf`",
    "`mardiaSkew`, `mardiaKurtosis`, invariance helpers",
    "`ggplot`, `aes`, `geom_line`, `labs`, `theme_bw`",
    "`cohens_d`",
    "`bootEGA`, `dimensionStability`",
    "`qgraph`, `EBICglasso`",
    "`NCT`, `plot`",
    "`estimateNetwork`, `bootnet`, `corStability`",
    "`mgm`, `predict`",
    "`boot.strength`, `averaged.network`, `bn.fit`, `predict`",
    "`confusionMatrix`",
    "`fct_reorder`",
    "`kable`",
    "`kable_styling`, `add_header_above`",
    "`cfa`, `fitMeasures`",
    "`semPaths`",
    "`corrplot`",
    "`read_sav`, `as_factor`",
    "`model_parameters`"
  )
)

# ---- Add Version and Loaded columns
packages_info <- pkgs %>%
  mutate(
    Version = vapply(Package, safe_version, character(1)),
    Loaded  = ifelse(vapply(Package, is_loaded, logical(1)), "Yes", "No")
  ) %>%
  select(Package, Version, Loaded, Functions)

# ---- Display table
knitr::kable(
  packages_info,
  caption = "Table 1. Packages used, versions, loaded status, and key functions"
) %>%
  kableExtra::kable_styling(full_width = FALSE, position = "center")
Table 1. Packages used, versions, loaded status, and key functions
Package Version Loaded Functions
readxl 1.4.3 No read_excel
tidyr 1.3.1 Yes pivot_longer, pivot_wider, reshaping support
dplyr 1.1.4 Yes select, mutate, filter, group_by, summarise, bind_rows
psych 2.4.1 Yes describe, polychoric, alpha, omega, splitHalf
semTools 0.5.6 Yes mardiaSkew, mardiaKurtosis, invariance helpers
ggplot2 3.5.2 Yes ggplot, aes, geom_line, labs, theme_bw
effectsize 0.8.9 Yes cohens_d
EGAnet 2.0.4 Yes bootEGA, dimensionStability
qgraph 1.9.8 Yes qgraph, EBICglasso
NetworkComparisonTest 2.2.2 Yes NCT, plot
bootnet 1.5.6 Yes estimateNetwork, bootnet, corStability
mgm 1.2.14 Yes mgm, predict
bnlearn 4.9 Yes boot.strength, averaged.network, bn.fit, predict
caret 6.0.94 Yes confusionMatrix
forcats 1.0.0 Yes fct_reorder
knitr 1.45 Yes kable
kableExtra 1.4.0 Yes kable_styling, add_header_above
lavaan 0.6.17 Yes cfa, fitMeasures
semPlot 1.1.6 No semPaths
corrplot 0.92 Yes corrplot
haven 2.5.4 Yes read_sav, as_factor
parameters 0.22.1 Yes model_parameters

Estamos utilizando 22 paquetes.

Design and generation of a formula for plotting centrality measures of two networks.

compareCentrality <- function(net1, net2,
                              include = c("Strength",
                                          "Closeness",
                                          "Betweenness",
                                          "ExpectedInfluence",
                                          "all",
                                          "All"),
                              orderBy = c("Strength",
                                          "Closeness",
                                          "Betweenness",
                                          "ExpectedInfluence"),
                              decreasing = T,
                              legendName = '',
                              net1Name = 'Network 1',
                              net2Name = 'Network 2'){
  
  if(include == "All" | include == "all"){
    include = c("Strength",
                "Closeness",
                "Betweenness",
                "ExpectedInfluence")
  }
  
  df <- centralityTable(net1, net2) %>% filter(measure %in% include)
  
  df %>% 
    mutate(graph = case_when(graph == 'graph 1' ~ net1Name,
                             graph == 'graph 2' ~ net2Name),
           graph = as.factor(graph),
           node = as.factor(node)) %>% 
    
    mutate(node = fct_reorder(node, value)) %>% 
    
    ggplot(aes(x = node, y = value, group = graph)) +
    
    geom_line(aes(linetype = graph), size = 1) +
    
    labs(x = '', y = '') +
    
    scale_linetype_discrete(name = legendName) +
    
    coord_flip() +
    
    facet_grid(~measure) +
    
    theme_bw()
  
}

1.2. Data Loading

The dataset from the #### project (Ethical Committee: ####) was loaded, including only the minimum set of variables required for the present study.

df_general <- read_sav("Estudio Telepresión - Fase 3 - ESO.sav")
df_general

1.3. Filtering

Filters were applied to remove participants who did not complete the isntrument or indicate their gender, and to exclude cases with missing values:

tele_df <- dplyr::select(df_general,c(Sexo,EdadV1,V2,AUTO:TPB))
names(tele_df)[1:3] <- c("sex","age","academic_year")
tele_df <- as.data.frame(tele_df)

tele_df$sex <- as.factor(case_when(tele_df$sex == 1 ~ "Male",
                         tele_df$sex == 2 ~ "Female",
                         tele_df$sex == 3 ~ "I do not know"))

tele_df$academic_year <- as.factor(case_when(tele_df$academic_year == 1 ~ "First academic period",
                         tele_df$academic_year == 2 ~ "First academic period",
                         tele_df$academic_year == 3 ~ "Second academic period",
                         tele_df$academic_year == 4 ~ "Second academic period"))

tele_df$PS <- as.factor(case_when(tele_df$PS == 0 ~ "Passive-Agressive",
                         tele_df$PS == 1 ~ "Assertive",
                         tele_df$PS == 2 ~ "Agressive",
                         tele_df$PS == 3 ~ "Passive"))

tele_df$TPB <- as.factor(case_when(tele_df$TPB == 1 ~ "Low risk",
                         tele_df$TPB == 2 ~ "Medium risk",
                         tele_df$TPB == 3 ~ "High risk"))

tele_df <- dplyr::select(tele_df, c("sex","academic_year","PS","TPB","age","AUTO","HETERO",
                                    "TP","SP","RS","CIUS","IDS","AB","VB","VCB","ACB","AF",
                                    "AA","APS"))

initial_rows <- nrow(tele_df)
tele_df <- na.omit(tele_df)

print(paste("Initial number of rows:", initial_rows, "| Final number of rows:", nrow(tele_df)))
## [1] "Initial number of rows: 820 | Final number of rows: 813"
tele_df

1.4. Variable Selection

Creation of dataframes tailored to meet the requirements of the forthcoming data analyses.

1. Dataframe for the analysis of individual items

The following table shows the first five rows of the dataset including the the variables used in the present study.

# Prepare the dataframe
network_variables <- dplyr::select(tele_df, age:APS)
pieColor <- rep("#EB9446", length(network_variables))

# Display first 5 rows with a formatted table
head(network_variables, 5) %>%
  kable(caption = "Table 2. Set of items from the WTM") %>%
  kable_styling(full_width = FALSE, position = "center")
Table 2. Set of items from the WTM
age AUTO HETERO TP SP RS CIUS IDS AB VB VCB ACB AF AA APS
14 10 27 3.500000 2.75 3.333333 2.533333 2.4 1.333333 1.333333 1.000000 1.000000 2.25 3.00 2.75
13 55 7 2.500000 2.25 2.000000 2.866667 3.2 1.666667 2.333333 1.000000 2.000000 1.00 1.25 1.25
13 55 53 1.833333 3.25 1.666667 1.333333 1.2 1.000000 2.000000 1.666667 1.000000 4.00 1.50 4.00
14 50 53 2.333333 2.75 2.000000 1.800000 1.0 1.666667 2.000000 1.333333 1.333333 4.00 4.00 3.75
14 45 47 2.500000 2.00 1.666667 2.600000 3.8 1.666667 1.333333 1.000000 1.000000 1.25 2.75 2.75

2. Dataframe for exploring sex differences

The following table shows the first five rows of the dataset including the the variables considered for the study and participants’ sex.

# Select CEI items, total score, and gender
NV_sex <- dplyr::select(tele_df, c(sex,age:APS))

# Display first 5 rows with styling
head(NV_sex, 5) %>%
  kable(caption = "Table 3. WTM items, total score, and gender variable") %>%
  kable_styling(full_width = FALSE, position = "center")
Table 3. WTM items, total score, and gender variable
sex age AUTO HETERO TP SP RS CIUS IDS AB VB VCB ACB AF AA APS
Female 14 10 27 3.500000 2.75 3.333333 2.533333 2.4 1.333333 1.333333 1.000000 1.000000 2.25 3.00 2.75
Female 13 55 7 2.500000 2.25 2.000000 2.866667 3.2 1.666667 2.333333 1.000000 2.000000 1.00 1.25 1.25
Female 13 55 53 1.833333 3.25 1.666667 1.333333 1.2 1.000000 2.000000 1.666667 1.000000 4.00 1.50 4.00
Male 14 50 53 2.333333 2.75 2.000000 1.800000 1.0 1.666667 2.000000 1.333333 1.333333 4.00 4.00 3.75
Female 14 45 47 2.500000 2.00 1.666667 2.600000 3.8 1.666667 1.333333 1.000000 1.000000 1.25 2.75 2.75

2. Dataframe for exploring academic year differences

The following table shows the first five rows of the dataset including the the variables considered for the study and participants’ academic year.

NV_academic_year <- dplyr::select(tele_df, c(academic_year,age:APS))

# Display first 5 rows with styling
head(NV_academic_year, 5) %>%
  kable(caption = "Table 3. WTM items, total score, and gender variable") %>%
  kable_styling(full_width = FALSE, position = "center")
Table 3. WTM items, total score, and gender variable
academic_year age AUTO HETERO TP SP RS CIUS IDS AB VB VCB ACB AF AA APS
First academic period 14 10 27 3.500000 2.75 3.333333 2.533333 2.4 1.333333 1.333333 1.000000 1.000000 2.25 3.00 2.75
First academic period 13 55 7 2.500000 2.25 2.000000 2.866667 3.2 1.666667 2.333333 1.000000 2.000000 1.00 1.25 1.25
First academic period 13 55 53 1.833333 3.25 1.666667 1.333333 1.2 1.000000 2.000000 1.666667 1.000000 4.00 1.50 4.00
First academic period 14 50 53 2.333333 2.75 2.000000 1.800000 1.0 1.666667 2.000000 1.333333 1.333333 4.00 4.00 3.75
First academic period 14 45 47 2.500000 2.00 1.666667 2.600000 3.8 1.666667 1.333333 1.000000 1.000000 1.25 2.75 2.75

2. DESCRIPTIVE ANALYSIS

2.1. Sociodemographic

Participants distribution by sex and age:

# sex distribution
sex_table <- table(tele_df$sex)
sex_percent <- round(prop.table(sex_table) * 100, 2)

sex_df <- data.frame(
  sex = names(sex_table),
  Frequency = round(as.numeric(sex_table), 2),
  Percentage = paste0(sex_percent, "%")
)

# Age summary
df_age <- tele_df[!is.na(tele_df$age), ]
age_summary <- summary(df_age$age)
age_sd <- round(sd(df_age$age), 2)

age_df <- data.frame(
  Statistic = names(age_summary),
  Value = round(as.numeric(age_summary), 2)
) %>%
  bind_rows(data.frame(Statistic = "Standard Deviation", Value = age_sd))

# Display tables
kable(sex_df, caption = "Table 5. sex distribution (frequency and %)") %>%
  kable_styling(full_width = FALSE, position = "center")
Table 5. sex distribution (frequency and %)
sex Frequency Percentage
Female 400 49.2%
I do not know 10 1.23%
Male 403 49.57%
kable(age_df, caption = "Table 6. Age summary statistics") %>%
  kable_styling(full_width = FALSE, position = "center")
Table 6. Age summary statistics
Statistic Value
Min. 12.00
1st Qu. 13.00
Median 14.00
Mean 13.84
3rd Qu. 15.00
Max. 18.00
Standard Deviation 1.24

2.2. Variables

Univariate descriptive analysis of the virables included in the analysis:

# Descriptive statistics
descrgen <- describe(network_variables)

# Build dataframe with selected stats
lgen <- list(c(1:15), descrgen$mean, descrgen$sd, descrgen$min,
             descrgen$max, descrgen$skew, descrgen$kurtosis)
lgen <- as.data.frame(lgen)

floor_pct  <- sapply(network_variables, function(x) mean(x == min(x, na.rm = TRUE), na.rm = TRUE) * 100)
ceiling_pct <- sapply(network_variables, function(x) mean(x == max(x, na.rm = TRUE), na.rm = TRUE) * 100)

lgen <- cbind(lgen,floor_pct)
lgen <- cbind(lgen,ceiling_pct)

# Round and rename columns
lgen <- lgen %>% 
  mutate_if(is.numeric, round, digits = 3)

colnames(lgen) <- c("Item", "Mean", "SD", "Min", "Max", "Skew", "Kurtosis", "Floor%", "Ceiling%")

# Display table
kable(lgen, caption = "Table 7. Descriptive statistics for WTM items") %>%
  kable_styling(full_width = FALSE, position = "center")
Table 7. Descriptive statistics for WTM items
Item Mean SD Min Max Skew Kurtosis Floor% Ceiling%
age 1 13.835 1.240 12 18.000 0.303 -0.415 15.744 0.246
AUTO 2 67.583 24.830 0 100.000 -0.574 -0.454 0.861 11.808
HETERO 3 59.513 26.603 0 100.000 -0.148 -0.871 1.353 11.070
TP 4 2.486 0.822 1 5.000 0.298 -0.229 4.182 0.246
SP 5 2.230 0.804 1 4.000 0.237 -0.770 10.578 2.583
RS 6 2.211 0.893 1 4.000 0.290 -0.984 15.867 4.797
CIUS 7 2.427 0.680 1 4.571 0.324 -0.101 0.984 0.123
IDS 8 1.490 0.825 1 5.000 2.105 4.044 52.399 0.492
AB 9 1.484 0.700 1 5.000 2.257 6.099 44.895 0.615
VB 10 1.860 0.899 1 5.000 1.362 1.614 24.477 1.353
VCB 11 1.394 0.705 1 5.000 2.642 8.209 59.410 0.984
ACB 12 1.280 0.595 1 5.000 3.380 13.940 67.651 0.492
AF 13 3.354 0.765 1 4.000 -1.316 1.046 2.091 36.039
AA 14 3.230 0.803 1 4.000 -0.762 -0.554 0.861 34.317
APS 15 3.290 0.745 1 4.000 -0.926 0.056 0.984 33.948

Assessment of multivariate normality assumptions based on Mardia’s tests for skewness and kurtosis, applied to the set of items included in the WTM instrument:

# Run Mardia's tests
skew_result <- mardiaSkew(network_variables)
kurt_result <- mardiaKurtosis(network_variables)

# Round and extract values
mardia_df <- data.frame(
  Test = c("Mardia Skewness", "Mardia Kurtosis"),
  Statistic = round(c(skew_result[2], kurt_result[2]), 3),
  p_value = round(c(skew_result[4], kurt_result[3]), 3)
)
row.names(mardia_df) <- 1:2
# Show table
kable(mardia_df, caption = "Table 8: Mardia's Tests for Multivariate Normality")
Table 8: Mardia’s Tests for Multivariate Normality
Test Statistic p_value
Mardia Skewness 5494.117 0
Mardia Kurtosis 48.142 0

2.3. Variables by sex

# Lista de ítems y etiquetas
items <- names(network_variables) #PARA MODIFICAR

item_labels <- c(
  "Age",
  "Autoassertivity",
  "Heteroassertivity",
  "Telepressure",
  "Supervision",
  "Restriction",
  "Problematic use of internet",
  "Suicide ideation",
  "Agression bullying",
  "Victim bullying",
  "Victim ciberbullying",
  "Agressor ciberbullying",
  "Family support",
  "Friends support",
  "Significant Person support"
) 

NV_sex <- NV_sex[NV_sex$sex != "I do not know",]

# Crear tabla de resultados
results <- lapply(seq_along(items), function(i) {
  item <- items[i]
  label <- item_labels[i]
  
  # Obtener estadísticas por grupo
  stats <- NV_sex %>%
    group_by(sex) %>%
    summarise(
      mean = round(mean(.data[[item]], na.rm = TRUE), 3),
      sd = round(sd(.data[[item]], na.rm = TRUE), 3)
    )
  
  male_val <- paste0(stats$mean[stats$sex == "Male"], " (", stats$sd[stats$sex == "Male"], ")")
  female_val <- paste0(stats$mean[stats$sex == "Female"], " (", stats$sd[stats$sex == "Female"], ")")
  
  # t-test
  t_result <- t.test(as.formula(paste(item, "~ sex")), data = NV_sex)
  t_stat <- round(t_result$statistic, 3)
  df <- round(t_result$parameter, 3)
  p_val <- round(t_result$p.value,3)
    
  # Efecto (Cohen's d)
  d_val <- cohens_d(as.formula(paste(item, "~ sex")), data = NV_sex)$Cohens_d
  d_val <- round(d_val, 3)
  
  # Combinar resultados
  data.frame(
    Item = label,
    Male = male_val,
    Female = female_val,
    T_statistic = paste0(t_stat, " (", df, ")"),
    p_value = p_val,
    Effect_size = d_val
  )
})

# Unir todos los ítems
final_table <- bind_rows(results)

# Mostrar tabla
kable(final_table, caption = "Table 9. Mean (SD) by sex, t-test Results, p-value and Effect Sizes for numeric variables")
Table 9. Mean (SD) by sex, t-test Results, p-value and Effect Sizes for numeric variables
Item Male Female T_statistic p_value Effect_size
Age 13.878 (1.251) 13.773 (1.214) -1.217 (800.604) 0.224 -0.086
Autoassertivity 71.998 (23.962) 63.163 (24.869) -5.126 (799.408) 0.000 -0.362
Heteroassertivity 65.089 (26.07) 53.8 (25.832) -6.164 (800.998) 0.000 -0.435
Telepressure 2.41 (0.827) 2.564 (0.807) 2.667 (800.759) 0.008 0.188
Supervision 2.181 (0.799) 2.294 (0.802) 2.008 (800.885) 0.045 0.142
Restriction 2.201 (0.91) 2.232 (0.877) 0.493 (800.316) 0.622 0.035
Problematic use of internet 2.364 (0.638) 2.48 (0.708) 2.437 (791.239) 0.015 0.172
Suicide ideation 1.363 (0.747) 1.597 (0.856) 4.138 (785.192) 0.000 0.292
Agression bullying 1.589 (0.747) 1.355 (0.581) -4.964 (757.796) 0.000 -0.350
Victim bullying 1.946 (0.959) 1.76 (0.804) -2.975 (779.264) 0.003 -0.210
Victim ciberbullying 1.403 (0.727) 1.367 (0.644) -0.729 (791.029) 0.466 -0.051
Agressor ciberbullying 1.319 (0.666) 1.228 (0.459) -2.264 (714.003) 0.024 -0.160
Family support 3.459 (0.679) 3.268 (0.819) -3.611 (772.53) 0.000 -0.255
Friends support 3.267 (0.776) 3.21 (0.812) -1.012 (798.818) 0.312 -0.071
Significant Person support 3.316 (0.727) 3.281 (0.74) -0.666 (800.519) 0.505 -0.047

2.4. Variables by academic year

# Lista de ítems y etiquetas
items <- names(network_variables) #PARA MODIFICAR

item_labels <- c(
  "Age",
  "Autoassertivity",
  "Heteroassertivity",
  "Telepressure",
  "Supervision",
  "Restriction",
  "Problematic use of internet",
  "Suicide ideation",
  "Agression bullying",
  "Victim bullying",
  "Victim ciberbullying",
  "Agressor ciberbullying",
  "Family support",
  "Friends support",
  "Significant Person support"
) 

# Crear tabla de resultados
results <- lapply(seq_along(items), function(i) {
  item <- items[i]
  label <- item_labels[i]
  
  # Obtener estadísticas por grupo
  stats <- NV_academic_year %>%
    group_by(academic_year) %>%
    summarise(
      mean = round(mean(.data[[item]], na.rm = TRUE), 3),
      sd = round(sd(.data[[item]], na.rm = TRUE), 3)
    )
  
  male_val <- paste0(stats$mean[stats$academic_year == "First academic period"], " (", stats$sd[stats$academic_year == "First academic period"], ")")
  female_val <- paste0(stats$mean[stats$academic_year == "Second academic period"], " (", stats$sd[stats$academic_year == "Second academic period"], ")")
  
  # t-test
  t_result <- t.test(as.formula(paste(item, "~ academic_year")), data = NV_academic_year)
  t_stat <- round(t_result$statistic, 3)
  df <- round(t_result$parameter, 3)
  p_val <- round(t_result$p.value,3)
    
  # Efecto (Cohen's d)
  d_val <- cohens_d(as.formula(paste(item, "~ academic_year")), data = NV_academic_year)$Cohens_d
  d_val <- round(d_val, 3)
  
  # Combinar resultados
  data.frame(
    Item = label,
    First_Period = male_val,
    Second_Period = female_val,
    T_statistic = paste0(t_stat, " (", df, ")"),
    p_value = p_val,
    Effect_size = d_val
  )
})

# Unir todos los ítems
final_table <- bind_rows(results)

# Mostrar tabla
kable(final_table, caption = "Table 10. Mean (SD) by academic year, t-test Results, p-value and Effect Sizes for numeric variables")
Table 10. Mean (SD) by academic year, t-test Results, p-value and Effect Sizes for numeric variables
Item First_Period Second_Period T_statistic p_value Effect_size
Age 13.03 (0.817) 14.911 (0.822) -32.377 (745.26) 0.000 -2.297
Autoassertivity 68.548 (24.376) 66.293 (25.403) 1.274 (730.568) 0.203 0.091
Heteroassertivity 60.686 (26.028) 57.945 (27.312) 1.444 (727.618) 0.149 0.103
Telepressure 2.434 (0.813) 2.556 (0.83) -2.101 (739.134) 0.036 -0.149
Supervision 2.439 (0.795) 1.949 (0.729) 9.12 (778.592) 0.000 0.638
Restriction 2.42 (0.874) 1.932 (0.842) 8.037 (761.665) 0.000 0.567
Problematic use of internet 2.392 (0.669) 2.472 (0.692) -1.66 (733.82) 0.097 -0.118
Suicide ideation 1.526 (0.854) 1.441 (0.782) 1.478 (778.999) 0.140 0.103
Agression bullying 1.453 (0.661) 1.526 (0.748) -1.435 (693.834) 0.152 -0.104
Victim bullying 1.883 (0.893) 1.83 (0.908) 0.834 (740.582) 0.405 0.059
Victim ciberbullying 1.351 (0.611) 1.452 (0.811) -1.957 (620.654) 0.051 -0.144
Agressor ciberbullying 1.212 (0.477) 1.371 (0.715) -3.591 (569.714) 0.000 -0.269
Family support 3.375 (0.761) 3.326 (0.77) 0.907 (742.657) 0.364 0.064
Friends support 3.169 (0.805) 3.31 (0.793) -2.492 (753.129) 0.013 -0.176
Significant Person support 3.282 (0.731) 3.301 (0.763) -0.363 (729.875) 0.717 -0.026

3. NETWORK ANALYSIS

3.1. Community Identification

A bootstrap exploratory graph analysis (EGA) was performed to estimate the number, structure, and stability of variable communities (i.e., network dimensions) based on their structural consistency. In this context, a community refers to a cluster of nodes that exhibit strong interconnections through edges.

3.1.1 Community Estimation

The analysis facilitates the identification and delineation of such communities, as well as the allocation of individual items to their corresponding dimensions.

set.seed(123)
communities <- bootEGA(network_variables,
                       model = "glasso",
                       type = "resampling",
                       iter = 500,
                       typicalStructure = TRUE,
                       plot.typicalStructure = TRUE)
Figure 1. Community Estimation

Figure 1. Community Estimation

3.1.2. Community Stability

Items were retained within their assigned communities only when their stability coefficients exceeded 0.70. Items falling below this threshold were excluded due to their potential to undermine the structural integrity of the corresponding dimension.

# Structural consistency
communities.dimstab <- dimensionStability(communities)
Figure 2. Community Stability

Figure 2. Community Stability

3.2. Gaussian Graphical Models (GGMs)

3.2.1. Network Estimation

# Definir los grupos (según índices de df_net)
group_list <- list(
  "Parenting Control" = c(1,5,6),
  "Social-Digital Use" = c(2,3,4,7),
  "Bullying" = c(9:12),
  "Social Support" = c(13:15),
  "Suicide"= 8
)

TP_glasso <- qgraph::qgraph(
  cor_auto(network_variables),
  graph        = "glasso",
  layout       = "spring",
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)
Figure 3. Gaussian Graphical Model Estimation

Figure 3. Gaussian Graphical Model Estimation

possible_edges <- (ncol(network_variables)*(ncol(network_variables)-1))/2
real_edges <- length(TP_glasso$Edgelist$weight)
density <- round(real_edges/possible_edges*100, digits = 2)

# Print explanation and result
cat("Number of potential edges:", possible_edges, "\n")
## Number of potential edges: 105
cat("Number of real edges:", real_edges, "\n")
## Number of real edges: 19
cat("Network density (percentage of actual edges out of all possible undirected edges):", density, "%\n")
## Network density (percentage of actual edges out of all possible undirected edges): 18.1 %
# Extract edge weights
edge_weights <- TP_glasso$Edgelist$weight

# Calculate summary statistics
mean_weight <- round(mean(edge_weights), digits = 3)
sd_weight <- round(sd(edge_weights), digits = 3)
max_weight <- round(max(edge_weights), digits = 3)
min_weight <- round(min(edge_weights), digits = 3)

# Print explanation and results clearly
cat("Standard Edge Weight \n")
## Standard Edge Weight
cat("The average edge weight was:", mean_weight, "(SD:", sd_weight, ")\n")
## The average edge weight was: 0.178 (SD: 0.298 )
cat("Minimum weight:", min_weight, "| Maximum weight:", max_weight, "\n")
## Minimum weight: -0.245 | Maximum weight: 0.574
# Calculate summary statistics
mean_weight <- round(mean(abs(edge_weights)), digits = 3)
sd_weight <- round(sd(abs(edge_weights)), digits = 3)
max_weight <- round(max(abs(edge_weights)), digits = 3)
min_weight <- round(min(abs(edge_weights)), digits = 3)

# Print explanation and results clearly
cat("\n")
cat("Absolute Values \n")
## Absolute Values
cat("The average edge weight was:", mean_weight, "(SD:", sd_weight, ")\n")
## The average edge weight was: 0.298 (SD: 0.169 )
cat("Minimum weight:", min_weight, "| Maximum weight:", max_weight, "\n")
## Minimum weight: 0.105 | Maximum weight: 0.574

3.2.2. Network characterization

TP_matrix <- as.matrix(network_variables)
invisible(capture.output({SV_mgm <-mgm(
    data = TP_matrix,
    type = rep("g", length(network_variables)),
    level = rep(1, length(network_variables)),
    k = 2,
    verbatim = TRUE,
    warnings = TRUE
  )
}))

#Predice la varianza explicada de cada nodo y otros parametros
SV_predic <- predict(SV_mgm, network_variables, error.continuous = "VarExpl")

TP_glasso_predicted <- qgraph::qgraph(
  cor_auto(network_variables),
  graph        = "glasso",
  layout       = "spring",
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = SV_predic$errors$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)
Figure 4. Gaussian Graphical Model Characterization

Figure 4. Gaussian Graphical Model Characterization

3.3. Centrality Measures

3.3.1. Centality Estimation

centrality_TP <- centralityPlot(TP_glasso, include = 
                                  c("Betweenness","Closeness","Strength","ExpectedInfluence"),
                                orderBy ="Betweenness", scale = "z-scores")
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
Figure 5. Centrality Measures Estimation

Figure 5. Centrality Measures Estimation

3.3.2. Centrality Stability

# 1. Estimar la red con EBICglasso
glasso_net <- estimateNetwork(network_variables, 
                              default = "EBICglasso", 
                              corMethod = "cor_auto",
                              verbose = FALSE)

# 2. Bootstrap de estabilidad (case-dropping)
set.seed(2025)
boot_results <- invisible(bootnet(
  estimateNetwork(network_variables, default = "EBICglasso", corMethod = "cor_auto"),
  type      = "case",
  nBoots    = 200,
  nCores    = 1,        # importante para ver mensajes
  statistics= c("strength", "closeness", "betweenness","expectedInfluence"),  # evita centralidades inestables
  verbose   = FALSE,
  maxErrors = 1000
))

# 3. Plot de estabilidad de centralidad (bootstrap case-dropping)
plot(boot_results,
     statistics = c("strength", "closeness", "betweenness","expectedInfluence"),
     labels = TRUE, order = "sample", legend = TRUE,
     color = "darkblue", line = TRUE,
     main = "Bootstrap Centrality Stability (Case Dropping)")
Figure 6. Centrality Measures Stability

Figure 6. Centrality Measures Stability

# 4. Coeficientes de estabilidad de centralidad
cs_coeffs <- corStability(boot_results, verbose = FALSE)

# Crear tabla
cs_table <- data.frame(
  Centrality = c("Betweenness", "Closeness", "Strength", "ExpectedInfluence"),
  CS_Coefficient = round(c(cs_coeffs[1], cs_coeffs[2], cs_coeffs[4], cs_coeffs[3]), 3)
)
row.names(cs_table) <- 1:4

# Mostrar como tabla en R Markdown
knitr::kable(cs_table, caption = "Table 10: Centrality Stability Coefficients (CS)",
align = c("r","r","r","r")) %>%
  kableExtra::kable_styling(full_width = FALSE, position = "center", bootstrap_options = c("striped","hover"))
Table 10: Centrality Stability Coefficients (CS)
Centrality CS_Coefficient
Betweenness 0.594
Closeness 0.594
Strength 0.750
ExpectedInfluence 0.750

4. NETWORK COMPARISON BY SEX

4.1. Network Comparison by age

TP_male <- NV_sex[NV_sex$sex == "Male",-1]
TP_female <- NV_sex[NV_sex$sex == "Female",-1]

set.seed(1234)
EN_male = estimateNetwork(TP_male, default = "EBICglasso", corMethod = "cor_auto")
EN_female = estimateNetwork(TP_female, default = "EBICglasso", corMethod = "cor_auto")
set.seed(1234)
invisible(capture.output({
  res <- NCT(
    TP_female,
    TP_male,
    p.adjust.methods = "BH",
    binary.data = FALSE,
    test.edges = TRUE,
    edges = "all",
    it = 1000
  )
}))

female_str <- round(res$glstrinv.sep[1], 3)
male_str <- round(res$glstrinv.sep[2], 3)
p_str <- ifelse(res$glstrinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$glstrinv.pval, 3)))

dif_net <- round(res$nwinv.real, 3)
p_net <- ifelse(res$nwinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$nwinv.pval, 3)))

p_edges <- paste(
  dplyr::filter(res$einv.pvals, `p-value` < 0.05) %>%
    dplyr::mutate(pair = paste(Var1, "-", Var2)) %>%
    dplyr::pull(pair),
  collapse = ", "
)

pvalue_edges <- paste(
  dplyr::filter(res$einv.pvals, `p-value` < 0.05)
)
pvalue_I4_I6 <- round(as.numeric(pvalue_edges[3]),3)
statistic_I4_I6 <- round(as.numeric(pvalue_edges[4]),3)

In terms of global strength, the female network showed a similary overall connectivity (S = 5.893) compared to the male network (S = 5.745). See Figure 7.

fmt_num <- function(x, d = 3) {
  if (is.numeric(x) && length(x) == 1 && is.finite(x)) sprintf(paste0("%.", d, "f"), x) else "NA"
}
fmt_p <- function(p, d = 3, eps = 1e-3) {
  if (is.numeric(p) && length(p) == 1 && is.finite(p)) format.pval(p, digits = d, eps = eps) else "NA"
}

plot(res, what="strength")
Figure 7. P-values for differences in network strength between sexs

Figure 7. P-values for differences in network strength between sexs

However, the network structure invariance test indicated a difference of 0.242 between the female and male networks. This difference was statistically significant (p = 0.049), indicating that the overall network structures differ (Figure 8). Follow-up edgewise tests, with p-values adjusted for multiple comparisons using the Benjamini–Hochberg procedure, showed that the behavior of the edge Item4–Item6 was statistically significant (E = NA, p = NA).

plot(res, what="network")
Figure 8. P-values for differences in network structure between sexs

Figure 8. P-values for differences in network structure between sexs

p_edges
## [1] ""

The Network Comparison Test did not identify any statistically significant differences in edge strength between the UN/NW and OWO networks after applying the Benjamini-Hochberg procedure to correct p-values for multiple comparisons. The significance level was set at the conventional threshold of 0.05.

4.2. Networks by sex

Female’s Network

# --- FEMALE GROUP ---

# Convert data to matrix
TP_female_matrix <- as.matrix(TP_female)

# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
  TP_mgm_female <- mgm(
    data = TP_female_matrix,
  type = rep("g", length(TP_female)),
  level = rep(1, length(TP_female)),
  k = 2,
  warnings = TRUE,
  verbatim = TRUE
  )
}))

# Predict explained variance (R²) and other node-level metrics
TP_predic_female <- predict(
  TP_mgm_female,
  TP_female,
  error.continuous = "VarExpl"
)

TP_glasso_predicted_female <- qgraph::qgraph(
  cor_auto(TP_female),
  graph        = "glasso",
  layout       = TP_glasso_predicted$layout,
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = SV_predic$errors$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)
Figure 9. Networks of Female Participants

Figure 9. Networks of Female Participants

Male’s network

# --- MALE GROUP ---

# Convert data to matrix
TP_male_matrix <- as.matrix(TP_male)

# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
  TP_mgm_male <- mgm(
    data = TP_male_matrix,
    type = rep("g", length(TP_male)),
    level = rep(1, length(TP_male)),
    k = 2,
    warnings = TRUE,
    verbatim = TRUE
  )
}))

# Predict explained variance (R²) and other node-level metrics
TP_predic_male <- predict(
  TP_mgm_male,
  TP_male,
  error.continuous = "VarExpl"
)

# Plot the network with R² represented as pie charts on nodes
TP_glasso_predicted_male <- qgraph::qgraph(
  cor_auto(TP_male),
  graph        = "glasso",
  layout       = TP_glasso_predicted$layout,
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = SV_predic$errors$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)
Figure 10. Networks of Male Participants

Figure 10. Networks of Male Participants

4.3. Centrality Measures

To evaluate the invariance of the centrality measures network across sexs.

# Compare centralities between networks
network_comparison <- compareCentrality(
  TP_glasso_predicted_female,
  TP_glasso_predicted_male,
  include = "all",
  legendName = "Centrality Measures by sex",
  net1Name = "Female",
  net2Name = "Male"
)
plot(network_comparison)
Figure 11. Centrality comparison between female and male networks

Figure 11. Centrality comparison between female and male networks

set.seed(1234)
invisible(capture.output({
  nct_test_centrality <- NCT(
    EN_female,
    EN_male,
    it = 1000,
    test.centrality = TRUE,
    p.adjust.methods = "BH",
    centrality = c("closeness", "betweenness", "strength", "expectedInfluence"),
  )
}))
# Extract p-values and convert to data frame
centrality_pvals <- as.data.frame(nct_test_centrality$diffcen.pval)

# Round to 3 decimals (and adjust "< 0.001" if you want a more elegant output)
centrality_pvals_rounded <- round(centrality_pvals, 3)

# Optional: add row names as a column
centrality_pvals_rounded$Node <- rownames(centrality_pvals_rounded)
centrality_pvals_rounded <- centrality_pvals_rounded[, c("Node", setdiff(names(centrality_pvals_rounded), "Node"))]
row.names(centrality_pvals_rounded) <- 1:15

# Display as table
kable(centrality_pvals_rounded, caption = "Table 11. P-values for Centrality Differences (NCT Test)") %>%
  kable_styling(full_width = FALSE, position = "center")
Table 11. P-values for Centrality Differences (NCT Test)
Node closeness betweenness strength expectedInfluence
age 0.673 1.000 0.885 0.885
AUTO 0.587 1.000 1.000 1.000
HETERO 0.587 1.000 0.742 1.000
TP 0.587 0.885 1.000 0.587
SP 0.687 1.000 0.587 0.587
RS 0.687 1.000 0.673 1.000
CIUS 0.587 0.587 0.060 0.885
IDS 0.687 1.000 0.673 0.673
AB 0.603 0.901 1.000 0.673
VB 0.741 1.000 0.776 0.587
VCB 0.741 1.000 0.645 0.587
ACB 0.687 1.000 0.587 0.741
AF 0.685 1.000 0.587 0.587
AA 0.776 1.000 0.687 1.000
APS 0.885 1.000 0.587 1.000

5. NETWORK COMPARISON BY ACADEMIC YEAR

5.1. Network Comparison by academic year

TP_first <- NV_academic_year[NV_academic_year$academic_year == "First academic period",-1]
TP_second <- NV_academic_year[NV_academic_year$academic_year == "Second academic period",-1]

set.seed(12)
EN_first = estimateNetwork(TP_first, default = "EBICglasso", corMethod = "cor_auto")
EN_second = estimateNetwork(TP_second, default = "EBICglasso", corMethod = "cor_auto")
set.seed(1234)
invisible(capture.output({
  res <- NCT(
    TP_second,
    TP_first,
    p.adjust.methods = "BH",
    binary.data = FALSE,
    test.edges = TRUE,
    edges = "all",
    it = 1000
  )
}))

second_str <- round(res$glstrinv.sep[1], 3)
first_str <- round(res$glstrinv.sep[2], 3)
p_str <- ifelse(res$glstrinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$glstrinv.pval, 3)))

dif_net <- round(res$nwinv.real, 3)
p_net <- ifelse(res$nwinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$nwinv.pval, 3)))

p_edges <- paste(
  dplyr::filter(res$einv.pvals, `p-value` < 0.05) %>%
    dplyr::mutate(pair = paste(Var1, "-", Var2)) %>%
    dplyr::pull(pair),
  collapse = ", "
)

pvalue_edges <- paste(
  dplyr::filter(res$einv.pvals, `p-value` < 0.05)
)
pvalue_I4_I6 <- round(as.numeric(pvalue_edges[3]),3)
statistic_I4_I6 <- round(as.numeric(pvalue_edges[4]),3)

In terms of global strength, the second network showed a similary overall connectivity (S = 5.243) compared to the first network (S = 5.779). See Figure 7.

fmt_num <- function(x, d = 3) {
  if (is.numeric(x) && length(x) == 1 && is.finite(x)) sprintf(paste0("%.", d, "f"), x) else "NA"
}
fmt_p <- function(p, d = 3, eps = 1e-3) {
  if (is.numeric(p) && length(p) == 1 && is.finite(p)) format.pval(p, digits = d, eps = eps) else "NA"
}

plot(res, what="strength")
Figure 7. P-values for differences in network strength between sexs

Figure 7. P-values for differences in network strength between sexs

However, the network structure invariance test indicated a statistically marginal difference (p = 0.512) of 0.146 between the second and first academic period networks (Figure 8). The posterior analysis of the edge do not detect any statistically significant difference to a local level among edges.

plot(res, what="network")
Figure 8. P-values for differences in network structure between sexs

Figure 8. P-values for differences in network structure between sexs

p_edges
## [1] ""

The Network Comparison Test did not identify any statistically significant differences in edge strength between the UN/NW and OWO networks after applying the Benjamini-Hochberg procedure to correct p-values for multiple comparisons. The significance level was set at the conventional threshold of 0.05.

4.2. Networks by Academic Period

second academic’s Network

# --- second GROUP ---

# Convert data to matrix
TP_second_matrix <- as.matrix(TP_second)

# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
  TP_mgm_second <- mgm(
    data = TP_second_matrix,
  type = rep("g", length(TP_second)),
  level = rep(1, length(TP_second)),
  k = 2,
  warnings = TRUE,
  verbatim = TRUE
  )
}))

# Predict explained variance (R²) and other node-level metrics
TP_predic_second <- predict(
  TP_mgm_second,
  TP_second,
  error.continuous = "VarExpl"
)

TP_glasso_predicted_second <- qgraph::qgraph(
  cor_auto(TP_second),
  graph        = "glasso",
  layout       = TP_glasso_predicted$layout,
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = SV_predic$errors$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)
Figure 9. Networks of second academic

Figure 9. Networks of second academic

first academic period’s network

# --- first GROUP ---

# Convert data to matrix
TP_first_matrix <- as.matrix(TP_first)

# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
  TP_mgm_first <- mgm(
    data = TP_first_matrix,
    type = rep("g", length(TP_first)),
    level = rep(1, length(TP_first)),
    k = 2,
    warnings = TRUE,
    verbatim = TRUE
  )
}))

# Predict explained variance (R²) and other node-level metrics
TP_predic_first <- predict(
  TP_mgm_first,
  TP_first,
  error.continuous = "VarExpl"
)

# Plot the network with R² represented as pie charts on nodes
TP_glasso_predicted_first <- qgraph::qgraph(
  cor_auto(TP_first),
  graph        = "glasso",
  layout       = TP_glasso_predicted$layout,
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = SV_predic$errors$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)
Figure 10. Networks of first academic period

Figure 10. Networks of first academic period

4.3. Centrality Measures

To evaluate the invariance of the centrality measures network across sexs.

# Compare centralities between networks
network_comparison <- compareCentrality(
  TP_glasso_predicted_second,
  TP_glasso_predicted_first,
  include = "all",
  legendName = "Centrality Measures by Academic Year",
  net1Name = "Second period",
  net2Name = "First period"
)
plot(network_comparison)
Figure 11. Centrality comparison between second and first networks

Figure 11. Centrality comparison between second and first networks

set.seed(1234)
invisible(capture.output({
  nct_test_centrality <- NCT(
    EN_second,
    EN_first,
    it = 1000,
    test.centrality = TRUE,
    p.adjust.methods = "BH",
    centrality = c("closeness", "betweenness", "strength", "expectedInfluence"),
  )
}))
# Extract p-values and convert to data frame
centrality_pvals <- as.data.frame(nct_test_centrality$diffcen.pval)

# Round to 3 decimals (and adjust "< 0.001" if you want a more elegant output)
centrality_pvals_rounded <- round(centrality_pvals, 3)

# Optional: add row names as a column
centrality_pvals_rounded$Node <- rownames(centrality_pvals_rounded)
centrality_pvals_rounded <- centrality_pvals_rounded[, c("Node", setdiff(names(centrality_pvals_rounded), "Node"))]
row.names(centrality_pvals_rounded) <- 1:15

# Display as table
kable(centrality_pvals_rounded, caption = "Table 11. P-values for Centrality Differences (NCT Test)") %>%
  kable_styling(full_width = FALSE, position = "center")
Table 11. P-values for Centrality Differences (NCT Test)
Node closeness betweenness strength expectedInfluence
age 0.455 1.000 0.325 1.000
AUTO 0.325 1.000 0.312 1.000
HETERO 0.325 1.000 1.000 1.000
TP 0.455 0.356 0.280 0.738
SP 0.455 0.710 1.000 1.000
RS 0.455 1.000 1.000 0.312
CIUS 0.280 0.455 0.280 1.000
IDS 0.321 0.781 0.455 0.587
AB 0.587 1.000 1.000 1.000
VB 0.671 1.000 0.455 1.000
VCB 1.000 1.000 1.000 0.582
ACB 1.000 0.240 0.587 0.653
AF 0.280 0.587 0.653 0.280
AA 0.312 1.000 0.672 1.000
APS 0.312 1.000 0.487 0.710

9. BAYESIAN NETWORK

The construction of the Bayesian network (BN) is carried out in two steps. Firstly, the estimation of the directed acyclic graph (DAG) is performed, and secondly, the BN model is fitted and validated with the study dataset.

5.1. DAG Estimation

For the DAG estimation phase, the PC Stable algorithm with no restrictions was employed. To ensure stability in the obtained DAG, a total of 200 bootstrap samples were drawn, and only the edges with a strength greater than 0.85 and a direction greater than 0.5 were retained in the final model.

set.seed(1234)
bootstr <- boot.strength(network_variables, R = 500, algorithm = "pc.stable")

avgnet <- averaged.network(bootstr, threshold = 0.85)

sp <- strength.plot(avgnet, bootstr, shape = "ellipse", render = FALSE)

TP_DAG_qgraph <- qgraph::qgraph(
  sp,
  layout       = "spring",
  groups       = group_list,
  #color        = group_colors,
  nodeNames    = item_labels,
  labels       = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)
Figure 12. Directed Acyclic Graph (DAG)

Figure 12. Directed Acyclic Graph (DAG)

A second model was generating including a blacklist to avoid IDS has influence in others variables.

# 1. Definir la blacklist: IDS no puede ser padre de nadie
vars <- names(network_variables)
destinos <- setdiff(vars, "IDS")

blacklist <- data.frame(
  from = "IDS",
  to   = destinos,
  stringsAsFactors = FALSE
)

# 2. Bootstrap con pc.stable usando la blacklist
set.seed(1234)
bootstr <- boot.strength(
  network_variables,
  R              = 500,
  algorithm      = "pc.stable",
  algorithm.args = list(
    blacklist = blacklist
  )
)

# 3. Red promedio con el umbral deseado
avgnet <- averaged.network(bootstr, threshold = 0.85)

# 4. Plot de fuerza
sp <- strength.plot(avgnet, bootstr, shape = "ellipse", render = FALSE)

TP_DAG_qgraph <- qgraph::qgraph(
  sp,
  layout       = "spring",
  groups       = group_list,
  #color        = group_colors,
  nodeNames    = item_labels,
  labels       = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  legend       = TRUE,
  legend.cex   = 0.8,
  layoutScale  = c(0.8, 0.8),
  rescale      = TRUE
)
Figure 12. Directed Acyclic Graph (DAG)

Figure 12. Directed Acyclic Graph (DAG)

5.2. Bayesian Network Validation

To fit and validate the BN model in the second phase of the process, the dataset was subdivided into 10 folds. A routine was implemented in which 90% of the folds were used to train the model, and the remaining 10% were used for testing. This cross-validation routine was repeated until all potential combinations were explore.

#Cross-validation (Check spatial cross-validation: blockCV)
set.seed(1234)
# Split data in 5 sets
kf <- dismo::kfold(nrow(network_variables), k = 10) # k-fold partitioning
models_fit <- c()

for (var in names(network_variables)) {
  kfold_fit <- c()
  aux_base <- data.frame('variable'=var)
  if(is.numeric(network_variables[[var]])){
    for(k in 1:10) {
      test <- network_variables[kf == k, ]
      train <- network_variables[kf != k, ]
      
      training <- bn.fit(avgnet,train)
      pred <- predict(training, node = var, data = test)
      
      z <- data.frame( R2 = R2(pred, test[[var]]),
                       RMSE = RMSE(pred, test[[var]]),
                       MAE = MAE(pred, test[[var]]))
      if(is.na(z$R2)){
        z[is.na(z$R2),] <- 0
      }
      kfold_fit <- rbind(kfold_fit, z) 
    }
    
    z <- apply(kfold_fit, 2, mean)
    z <- cbind(aux_base,as.data.frame(t(z)))
    models_fit <- rbind(models_fit, z)
  }else{
    for(k in 1:10) {
      # Split data in test and train
      test <- network_variables[kf == k, ]
      train <- network_variables[kf != k, ]
      
      training <- bn.fit(avgnet,train)
      pred <- predict(training, node = var, data = test)
      
      validation <- confusionMatrix(as.factor(pred),test[[var]])
      z1 <- as.data.frame(t(validation$overall))
      
      if(length(levels(network_variables[[var]])) > 2){
        z2 <- as.data.frame(t(sapply(as.data.frame(validation$byClass),mean))) 
      }else{
        z2 <- as.data.frame(t(sapply(validation$byClass,mean))) 
      }
      z <- cbind(z1, z2)
      
      kfold_fit <- rbind(kfold_fit, z)
    }
    z <- apply(kfold_fit, 2, mean)
    z <- cbind(aux_base,z)
    models_fit <- rbind(models_fit, z)
  }
}

# Plot the network with R² represented as pie charts on nodes
TP_BN_qgraph <- qgraph::qgraph(
  sp,
  layout       = "spring",
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = models_fit$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)
Figure 13. Bayesian Network (BN) Model

Figure 13. Bayesian Network (BN) Model

10. PREDICTAIBILITY CAPACITY

# GGM - General
ggm_general <- SV_predic$errors %>%
  transmute(
    Variable = items,
    R2_gen = round(R2, 3),
    RMSE_gen = round(RMSE, 3)
  )

# GGM - Female
ggm_female <- TP_predic_female$errors %>%
  transmute(
    R2_fem = round(R2, 3),
    RMSE_fem = round(RMSE, 3)
  )

# GGM - Male
ggm_male <- TP_predic_male$errors %>%
  transmute(
    R2_male = round(R2, 3),
    RMSE_male = round(RMSE, 3)
  )

# GGM - First
ggm_first <- TP_predic_first$errors %>%
  transmute(
    R2_first = round(R2, 3),
    RMSE_first = round(RMSE, 3)
  )

# GGM - Second
ggm_second <- TP_predic_second$errors %>%
  transmute(
    R2_second = round(R2, 3),
    RMSE_second = round(RMSE, 3)
  )

# BN model
bn <- models_fit %>%
  transmute(
    R2_bn = ifelse(R2 == 0, "-", round(R2, 3)),
    RMSE_bn = ifelse(RMSE == 0, "-", round(RMSE, 3))
  )

# Combine all into final table
final_table <- bind_cols(
  ggm_general, 
  ggm_female, 
  ggm_male, 
  ggm_first, 
  ggm_second,
  bn
)

# Set column names
names(final_table) <- c(
  "Variable",
  rep(c("R²", "RMSE"), times = 6)  # General, Female, Male, First, Second, BN
)

# Generate final table
final_table %>%
  kable(align = "lcccccccccccc", booktabs = TRUE,
        caption = "Table 12. Explained variance and RMSE of partial-correlation network (GGM) and BN models") %>%
  add_header_above(c(" " = 1, 
                     "General" = 2, 
                     "Female" = 2, 
                     "Male" = 2,
                     "First Period" = 2,
                     "Second Period" = 2,
                     "DAG" = 2)) %>%
  add_header_above(c(" " = 1, 
                     "GGMs" = 10, 
                     "BN" = 2)) %>%
  kable_styling(full_width = FALSE, position = "center")
Table 12. Explained variance and RMSE of partial-correlation network (GGM) and BN models
GGMs
BN
General
Female
Male
First Period
Second Period
DAG
Variable RMSE RMSE RMSE RMSE RMSE RMSE
age 0.166 0.912 0.180 0.905 0.146 0.923 0.102 0.946 0.057 0.970 0.116 1.178
AUTO 0.556 0.666 0.533 0.683 0.552 0.668 0.579 0.648 0.536 0.680 0.457 18.265
HETERO 0.483 0.718 0.433 0.752 0.507 0.701 0.473 0.725 0.488 0.715
TP 0.298 0.837 0.382 0.785 0.248 0.866 0.244 0.868 0.385 0.783
SP 0.452 0.740 0.436 0.750 0.461 0.733 0.371 0.792 0.460 0.734 0.415 0.617
RS 0.392 0.779 0.351 0.804 0.431 0.753 0.312 0.829 0.442 0.746
CIUS 0.431 0.754 0.536 0.681 0.335 0.814 0.458 0.736 0.434 0.751 0.38 0.538
IDS 0.373 0.791 0.414 0.765 0.323 0.822 0.412 0.766 0.336 0.814 0.359 0.659
AB 0.555 0.666 0.563 0.660 0.554 0.667 0.508 0.701 0.619 0.617
VB 0.528 0.687 0.601 0.631 0.507 0.701 0.509 0.700 0.561 0.661 0.476 0.652
VCB 0.511 0.699 0.412 0.766 0.578 0.648 0.438 0.749 0.586 0.643
ACB 0.573 0.653 0.451 0.740 0.637 0.602 0.465 0.730 0.666 0.577 0.518 0.406
AF 0.410 0.768 0.353 0.804 0.516 0.695 0.410 0.767 0.454 0.738 0.29 0.642
AA 0.447 0.743 0.512 0.697 0.390 0.780 0.463 0.732 0.408 0.768 0.431 0.607
APS 0.527 0.688 0.545 0.674 0.527 0.687 0.537 0.679 0.544 0.674